//options
clear all
set maxvar  32000, permanently
set matsize 11000, permanently
set more off, permanently
global bootstraps 1000
set seed 23

// macros
global klmshare				: env klmshare
global projects				: env projects
global klmperry                         : env klmperry

global perrydatas       = "${klmperry}/CBA/data/perry/raw"
global perrydata        = "${klmperry}/CBA/data/perry/clean"
global masterinter  	= "${klmshare}/CurrentRAs/FredB/test"
global output           = "${projects}/ece_parenting/tex_files/other"
global dataperry        = "$klmperry/PerryPreschool/Data/Perry_PARI_and_Other_Data/PARI/DATA_PARI/"
global abcjpeanalysis   = "${klmshare}/Data_Central/Abecedarian/data/ABC-CARE/extensions/cba-iv/"
global dataihdp         = "${projects}/ece_parenting/data"

// base
cd  $perrydatas
use perry_base.dta, clear
gen     treat = 1 if C2 == 2
replace treat = 0 if C2 == 1
gen perry = 1
keep    id treat sb_iq_5 perry
rename sb_iq_5 earlyiq
tempfile perry
save "`perry'", replace

// abc
cd $abcjpeanalysis
use append-abccare_iv.dta, clear
keep if program == "abc"
egen sb45y = rowmean(sb4y sb5y)
rename sb45y earlyiq
keep id treat earlyiq
gen abc = 1
tempfile abc
save "`abc'", replace

// ihdp
cd  $dataihdp
use ihdp_data, clear
rename _all, lower

replace   bw = bw/1000
destring bwg , replace 

// impute iq
replace   stndscor =  iqcage     if stndscor ==.
replace   iqcage   =  stndscor   if iqcage ==.
foreach var of varlist stndscor iqcage { 
	replace `var' = mdi_24cor    if stndscor ==. & iqcage ==.
}
replace   stndscor =  iqcage     if stndscor ==.
replace   iqcage   =  stndscor   if iqcage ==.

// list relevant variables
global baseline_child      twin sex bw anga black hispanic 
global baseline_mother     mage meduc works married 
global baseline_household  welfare tot_siblings_natural employed_adult
global baseline_economy    employment medinc gpc
global outputs             iqcage stndscor 
global inputs              cum_avg_daycare_36m_sum

// mark baseline sample
reg  $baseline_child $baseline_mother $baseline_household $baseline_economy $outputs $inputs
gen  sample1 = e(sample) 
keep if sample1 == 1

// parenting latents, ages 1 and 3
foreach var of varlist alt_subscale_1_12-alt_subscale_6_12 alt_subscale_1_36-alt_subscale_8_36 {
	summ `var'
	gen  `var'_std = (`var' - r(mean))/r(sd)
}
# delimit
sem (_cons@0 X -> alt_subscale_1_12_std) (_cons@0 X -> alt_subscale_2_12_std) (_cons@0 X -> alt_subscale_3_12_std) 
	(_cons@0 X -> alt_subscale_4_12_std) (_cons@0 X -> alt_subscale_5_12_std) (_cons@0 X -> alt_subscale_6_12_std), method(adf);
predict parenting_age1, latent;
sem (_cons@0 X -> alt_subscale_1_36_std) (_cons@0 X -> alt_subscale_2_36_std) (_cons@0 X -> alt_subscale_3_36_std) (_cons@0 X -> alt_subscale_4_36_std)
	(_cons@0 X -> alt_subscale_5_36_std) (_cons@0 X -> alt_subscale_6_36_std) (_cons@0 X -> alt_subscale_7_36_std) (_cons@0 X -> alt_subscale_8_36_std), 
																																						 cov(e.alt_subscale_1_36_std*e.alt_subscale_2_36_std)
																																						 cov(e.alt_subscale_3_36_std*e.alt_subscale_4_36_std)
																																						 cov(e.alt_subscale_5_36_std*e.alt_subscale_6_36_std)
																																						 cov(e.alt_subscale_7_36_std*e.alt_subscale_8_36_std)
																																						 cov(e.alt_subscale_4_36_std*e.alt_subscale_7_36_std)
																																						 cov(e.alt_subscale_2_36_std*e.alt_subscale_7_36_std) method(adf);		
 predict parenting_age3, latent;
# delimit cr
egen    parenting_ages13 = rowmean(parenting_age1 parenting_age3)

// standardize latents in sample
foreach var of varlist parenting_ages13 {
	summ    `var' if tg == 0
	gen     `var'_std = (`var' - r(mean))/r(sd)
	replace `var'_std = `var'_std + 100
}
// log and standardize childcare variable 
summ    cum_avg_daycare_36m_sum  if tg == 0 
replace cum_avg_daycare_36m_sum  = (cum_avg_daycare_36m_sum    - r(mean))/r(sd)
replace cum_avg_daycare_36m_sum  =  cum_avg_daycare_36m_sum    + r(mean)

global inputs_std cum_avg_daycare_36m_sum parenting_ages13_std
keep   ihdp site tg bwg $baseline_child $baseline_mother $baseline_household $baseline_economy $inputs_std $outputs
rename ihdp id
rename tg treat
rename iqcage earlyiq
gen    ihdp = 1 
keep   id treat earlyiq ihdp twin bwg site 
append using "`perry'"
append using "`abc'"

// specs for run
gen ihdps = 1 if ihdp == 1 & twin == 0
gen ihdpt = 1 if ihdp == 1 & twin == 1
gen ihdpp = 1 if ihdp == 1 & twin != .

// strata
global perry_strata
global abc_strata
global ihdps_strata strata(bwg site treat)
global ihdpt_strata strata(bwg site treat)
global ihdpp_strata strata(bwg site treat)

// treatment effecs
foreach prog in abc perry ihdps ihdpt ihdpp {
	bootstrap, ${`prog'_strata} reps($bootstraps): reg earlyiq treat if `prog' == 1
	matrix b = e(b)
	matrix V = e(V)
	matrix se`prog' = sqrt(V[1,1])
	matrix md`prog' = b[1,1]
	matrix  t`prog' = abs(md`prog'[1,1]/se`prog'[1,1])
	matrix df`prog' = e(N) - e(rank)
	matrix  p`prog' = 2*(1 - normal(t`prog'[1,1]))
	matrix  n`prog' = e(N)
	matrix `prog' = [md`prog',p`prog',n`prog']
}

foreach prog in abc perry ihdps ihdpt ihdpp {
	summ    earlyiq                     if treat == 0 & `prog' == 1
	replace earlyiq = earlyiq - r(mean) if `prog' == 1
}

// two samples ts
foreach var in s t p {
	replace ihdp`var' = 0 if abc == 1 | perry == 1
}

foreach prog in abc perry {
	foreach sib in s t p {
		bootstrap, reps($bootstraps): reg earlyiq ihdp`sib' if (`prog' == 1 | ihdp`sib' == 1) & treat == 1
		matrix b = e(b)
		matrix V = e(V)
		matrix se`prog'_s`sib' = sqrt(V[1,1])
		matrix md`prog'_s`sib' = b[1,1]
		matrix  t`prog'_s`sib' = abs(md`prog'_s`sib'[1,1]/se`prog'_s`sib'[1,1])
		matrix df`prog'_s`sib' = e(N) - e(rank)
		matrix  p`prog'_s`sib' = 2*(1 - normal(t`prog'_s`sib'[1,1]))
	}
	matrix `prog' = [`prog',p`prog'_sp,p`prog'_ss,p`prog'_st]
	matrix colnames `prog' = md p n pip pis pit
}

// put together and plot
foreach sib in s t p {
	matrix ihdp`sib' = [ihdp`sib',.,.,.]
}

matrix all = J(1,6,.)
foreach prog in ihdpp ihdps ihdpt perry abc {
	matrix all = [all \ `prog']
}
matrix all = all[2...,1...]
matrix colnames all = md p n pip pis pit

clear
svmat all, names(col)

// mean difference
tostring  md, gen(mdstr) force format(%15.1fc) 
replace   mdstr = "{&Delta}{sup:ABC}   = " + mdstr if _n == 5
replace   mdstr = "{&Delta}{sup:Perry} = " + mdstr if _n == 4
replace   mdstr = "{&Delta}{sup:IHDP} = "  + mdstr if _n == 1
replace   mdstr = "{&Delta}{sup:IHDP} = "  + mdstr if _n == 2
replace   mdstr = "{&Delta}{sup:IHDP} = "  + mdstr if _n == 3

gen ptp1 = 19.5
gen ptp2 = 18.5

gen ptp3 = 16
gen ptp4 = 15

gen ptp5 = 12.5
gen ptp6 = 11.5


foreach sib in s t p {
	gen      ho`sib' = pi`sib'
	tostring ho`sib', replace force format(%15.2fc)
	replace  ho`sib' = "{it:p} = " + ho`sib'
}
gen      hoperry = "H{sub:0}: {&Delta}{sup:Perry} = {&Delta}{sup:IHDP}"
gen      hoabc   = "H{sub:0}: {&Delta}{sup:ABC} = {&Delta}{sup:IHDP}"

// n
tostring n, replace force format(%15.0fc)
replace  n = "N = " + n

foreach num of numlist 1(1)10 {
	gen n`num' = `num'
}


#delimit
twoway     (bar     md    n2    if _n == 4, color(gs10) barwidth(.9)       yline(0, lpattern(solid) lwidth(vthick) lcolor(gs14)))
	   (scatter ptp1  n2    if _n == 4, msymbol(none) mlabel(mdstr)    mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   (scatter ptp2  n2    if _n == 4, msymbol(none) mlabel(n)        mlabposition(12) mlabcolor(gs0) mlabsize(small))  
	   (scatter md    n2    if _n == 4 & p    < 0.05, msize(medlarge)  mcolor(gs0) msymbol(circle))
	   (scatter md    n2    if _n == 4 & p >= .05 & p    <= 0.10, msize(medlarge)  mcolor(gs0) msymbol(square))
	  
	   (bar     md    n3    if _n == 5, color(gs10) barwidth(.9)       yline(21, lpattern(solid) lcolor(gs14)))
	   (scatter ptp1  n3    if _n == 5, msymbol(none) mlabel(mdstr)    mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   (scatter ptp2  n3    if _n == 5, msymbol(none) mlabel(n)        mlabposition(12) mlabcolor(gs0) mlabsize(small))  
	   (scatter md    n3    if _n == 5 & p    < 0.05, msize(medlarge)  mcolor(gs0) msymbol(circle))
	   (scatter md    n3    if _n == 5 & p >= .05 & p    <= 0.10, msize(medlarge)  mcolor(gs0) msymbol(square))
	   
	   (scatter ptp3  n5 if _n == 4, msymbol(none) mlabel(hoperry)  mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   (scatter ptp4  n5 if _n == 4, msymbol(none) mlabel(hop)      mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   
	   (scatter ptp5  n5 if _n == 5, msymbol(none) mlabel(hoabc)    mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   (scatter ptp6  n5 if _n == 5, msymbol(none) mlabel(hop)      mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   
	   (bar     md    n5    if _n == 1, color(gs10) barwidth(.9))
	   (scatter ptp1  n5    if _n == 1, msymbol(none) mlabel(mdstr)    mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   (scatter ptp2  n5    if _n == 1, msymbol(none) mlabel(n)        mlabposition(12) mlabcolor(gs0) mlabsize(small))  
	   (scatter md    n5    if _n == 1 & p    < 0.05, msize(medlarge)  mcolor(gs0) msymbol(circle))
	   (scatter md    n5    if _n == 1 & p >= .05 & p    <= 0.10, msize(medlarge)  mcolor(gs0) msymbol(square))
	   
	   (scatter ptp3  n6 if _n == 4, msymbol(none) mlabel(hoperry)  mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   (scatter ptp4  n6 if _n == 4, msymbol(none) mlabel(hos)      mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   
	   (scatter ptp5  n6 if _n == 5, msymbol(none) mlabel(hoabc)    mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   (scatter ptp6  n6 if _n == 5, msymbol(none) mlabel(hos)      mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   
	   (bar     md    n6    if _n == 2, color(gs10) barwidth(.9))
	   (scatter ptp1  n6    if _n == 2, msymbol(none) mlabel(mdstr)    mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   (scatter ptp2  n6    if _n == 2, msymbol(none) mlabel(n)        mlabposition(12) mlabcolor(gs0) mlabsize(small))  
	   (scatter md    n6    if _n == 2 & p    < 0.05, msize(medlarge)  mcolor(gs0) msymbol(circle))
	   (scatter md    n6    if _n == 2 & p >= .05 & p    <= 0.10, msize(medlarge)  mcolor(gs0) msymbol(square))
	   
	   (scatter ptp3  n7 if _n == 4, msymbol(none) mlabel(hoperry)  mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   (scatter ptp4  n7 if _n == 4, msymbol(none) mlabel(hot)      mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   
	   (scatter ptp5  n7 if _n == 5, msymbol(none) mlabel(hoabc)    mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   (scatter ptp6  n7 if _n == 5, msymbol(none) mlabel(hot)      mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   
	   (bar     md    n7    if _n == 3, color(gs10) barwidth(.9))
	   (scatter ptp1  n7    if _n == 3, msymbol(none) mlabel(mdstr)    mlabposition(12) mlabcolor(gs0) mlabsize(small))
	   (scatter ptp2  n7    if _n == 3, msymbol(none) mlabel(n)        mlabposition(12) mlabcolor(gs0) mlabsize(small))  
	   (scatter md    n7    if _n == 3 & p    < 0.05, msize(medlarge)  mcolor(gs0) msymbol(circle))
	   (scatter md    n7    if _n == 3 & p >= .05 & p    <= 0.10, msize(medlarge)  mcolor(gs0) msymbol(square))
	   
	   ,
		  legend(label(1 "{&Delta}: Treatment Effect") label(4 "{&Delta} {&ne} 0 ({it:p}-value < 0.05)")
			 label(5 "{&Delta} {&ne} 0 ({it:p}-value {&isin} [0.05,.10])")
		        cols(3) rows(1) order(4 5) size(medsmall) position(6) region(lcolor(gs0)))
				 
		  ylabel(0[7]21, labsize(medium) angle(h) glcolor(gs14) glpattern(solid))  
		  xlabel(2 "{bf:Perry}" 3 "{bf:ABC}" 5 "Pooled" 6 `" "Singletons" "{bf:IHDP}" "' 7 "Twins",  labsize(medlarge) grid glcolor(gs14) glpattern(solid))
		  ytitle("{&Delta}: Treatment Effect", size(medlarge))
		  xtitle("")
		  graphregion(color(white)) plotregion(fcolor(white)); 
#delimit cr
cd $output
graph export cskills_alltwinning.eps, replace
